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ABSTRAGT 

The formation history of rich clnsters is investigated using a hybrid N-body simulation 
in which high spatial and mass resolution can be achieved self-consistently within a small 
region of a very large volume. The evolution of three massive clusters is studied via mass 
accretion, spherically-averaged density profiles, three-dimensional and projected shapes, 
and degree of substructure. Each cluster consists of ~ 4 x 10® particles at the present 
epoch and in the case that rich cluster evolution is well-described by a 1-parameter family, 
the simulations have sufficient resolution to demonstrate this. At z — 0 the clusters have 
similar masses, M(r < 1.5/i“^Mpc) ~ 2 x 10A^h~^MQ, and similar spherically-averaged 
density profiles, however markedly different formation histories are observed. No single, 
dominant pattern is apparent in the time variation of the mass accretion rate, the cluster 
shape, or the degree of substructure. Although not a statistically large sample, these 
objects suggest that the detailed formation history of rich clusters cannot be characterized 
by a simple 1-parameter family. These results suggest that the use of observations of rich 
clusters over a wide range of redshifts to constrain cosmological parameters may not be 
entirely straightforward. 

Subject headings: cosmology: dark matter — cosmology: large-scale structure of the uni¬ 
verse — cosmology: theory — methods: numerical 
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1. INTRODUCTION 


The evolution history of clusters of galaxies is a potentially powerful constraint on 
theories of the formation of large-scale structure in the universe. Although a statistically 
large sample does not yet exist, it is possible to identify clusters to high redshift (eg. Smail 
& Dickinson 1995; Bower & Smail 1997; Deltorn et al. 1997; Luppino & Kaiser 1997; Steidel 
et al. 1997) and it is hoped that eventually a complete description of the time evolution 
of these objects will be obtained. Additionally, detailed studies of clusters have yielded 
evidence for signihcant amounts of substructure in many clusters, even those which appear 
smooth and round in projection (eg. Beers & Geller, 1983; Jones & Forman 1984; Dressier 
& Shectman 1988; West & Bothun 1990; Davis & Mushotzky 1993; Miyaji et al. 1993; 
Mushotzky 1993; White, Briel & Henry 1993; Bird 1994a,b; Zabludoff & Zaritsky 1995). 
Optical, X-ray, and kinematic evidence has been found for recent mergers of a number of 
clusters with smaller systems. It appears that at least one third of all clusters are not fully 
relaxed systems and it is possible that many clusters are still in the process of assembling 
even today. The fraction of clusters containing signihcant amounts of substructure at the 
present epoch is a potentially powerful constraint on the value of the density parameter 
(eg. Richstone, Loeb & Turner 1992; Bartelmann, Ehlers & Schneider 1993; Kauffmann 
& White 1993; Lacey & Cole 1993; Mohr et al. 1995) and, therefore, cluster substructure 
investigations are of considerable interest. 

It has been hoped that numerical simulations of cluster formation would provide useful 
constraints on large-scale structure theories via comparisons of the formation histories of 
simulated clusters and the observed cluster population. Direct comparisons are, however, 
problematical for a number of reasons. First, pure N-body simulations follow only the 
evolution of the dominant, dissipationless mass component of the universe, neglecting hy¬ 
drodynamics. In such simulations a direct comparison of a theoretical distribution of light 
(eg. galaxies and X-ray gas) to that of observed clusters is not possible. Under the assump¬ 
tion that density peaks of an appropriate mass scale correspond to the likely sites of galaxy 
formation, however, it is possible to locate groups of particles within the simulations that 
may be associated fairly with the dark matter halos of individual galaxies. Additionally, 
from studies of the coherent weak shear held associated with gravitational lens clusters, 
the mass of the clusters is certainly dominated by dark matter and in addition it appears 
that the mass distribution within the clusters traces the smoothed light distribution quite 
well (eg. Bonnet, Mellier & Fort 1994; Fahlman et al. 1994; Smail et al. 1995; Kneib et 
al. 1996; Seitz et al. 1996; Squires et al. 1996ab; Smail et al. 1997). Therefore, cluster 
simulations which include only a dark matter component should yield fairly reasonable 
results for comparison with observation, at least to a good hrst approximation. 

The worst problem to plague simulations of cluster formation is simply one of reso¬ 
lution, both in terms of the gravitational force calculation on small scales as well as the 
mass per particle. That is, within the cluster enviro nm ent itself, one would like to resolve 
the physical scales associated with galaxies (distances of order a few tens of kiloparsecs, 
using a mass per particle of order IO^Mq). Ideally, of course, one would like to achieve 
such resolution inside a simulation volume which is itself a “fair sample” of the universe 
(of order 10^/i“^Mpc^). Current computing platforms, however, do not allow such a high 
level of resolution throughout a large simulation volume. 
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Therefore, high-resolution simulations of cluster formation often follow a scenario in 
which a simulation of a large volume of the universe is hrst run at moderate resolution. 
The hnal timestep of this simulation is then searched for peaks in the smoothed mass 
density which would correspond fairly to clusters of galaxies. A sphere of a given radius 
(typically r ~ 10h~^ Mpc) centered on the peak is then excised from the initial conditions 
of the simulation and populated with a large number of particles of small (sub-galactic) 
mass (eg. Bromley et ah, 1995; Carlberg, 1994). These smaller peak particles are then 
evolved from the initial timestep to the present epoch, subject to an external potential 
which is intended to model the correct tidal held due to the local universe that the cluster 
would experience as it evolves. DifRculties with such simulations are insuring: (1) the 
radius of the sphere is large enough to include all the mass that should be accreted by the 
cluster up to the present epoch and (2) the model external potential fairly represents the 
actual tidal held the cluster would experience if one simply ran the entire simulation at an 
unachieveably high resolution level. 

Here we simulate the formation of 3 rich clusters at very high resolution and inves¬ 
tigate the similarities and diherences of their evolution histories. All three clusters have 
similar masses at the present day and for simplicity a standard cold dark matter universe 
is adopted. The clusters form inside a large computational volume (8.0 x lO^h ^Mpc^) and 
high resolution is achieved without the use of either constrained initial conditions or the 
excision of peaks from a large-scale density held. Instead, a hybrid N-body code utilized. 
This particular code allows high spatial and mass resolution to be achieved simultaneously 
within small selected regions of a very large primary simulation volume. High resolution 
is obtained by nesting small simulations self-consistently inside larger simulations, result¬ 
ing in a “power zoom” ehect. The basic premise behind the N-body code used for the 
investigation is outlined in §2. Details of the simulations performed are summarized in §3, 
results of the analysis of the simulations are presented in §4, and a discussion of the results 
is given in §5. 

2. HIERARCHICAL PARTICLE MESH (HPM) 

The code used to run the simulations is the Hierarchical Particle Mesh (HPM) code 
written by J. V. Villumsen (Villumsen, 1989). The heart of the HPM code is a standard 
particle mesh (PM) cosmological simulation in which mass density is assigned to a grid 
using a cloud-in-cell (CIC) weighting scheme and Poisson’s equation is solved using fast 
Fourier techniques. Although very fast, PM simulations suffer from limited spatial resolu¬ 
tion, the force being softer than Newtonian on scales smaller than about 2 grid cells. In 
order to gain both spatial and mass resolution in a small region of the primary simulation 
volume, the HPM code allows small simulations (“subgrids”) to be nested self-consistently 
within the main simulation. By nesting subgrids inside subgrids one can progressively 
build up to very high resolution in a limited region of the total simulation volume. Details 
of the force calculations and the generation of initial conditions for multi-grid simulations 
are given in Villumsen (1989); here we present only a brief outline of the premise behind 
HPM. 

It is important to note that a multi-grid simulation using HPM is an iterative process. 
To begin, an ordinary PM simulation of a large volume of the universe is run from the 


3 



desired starting epoch {z ~ 30) to the present epoch. Periodic bonndary conditions are 
imposed on this grid and the simnlation is carried ont in a manner similar to all conven¬ 
tional PM simnlations of the formation of large-scale strnctnre. Thronghont, this large 
grid shall be referred to as the “top grid”; it constitntes the fnndamental compntational 
volnme of the simnlation. 

A small region of interest which is to be rnn in high resolntion mode (eg. a clnster 
environment) is then selected from the final timestep of the top grid. Using the previonsly 
recorded timesteps of the top grid calcnlation, those top grid particles which pass throngh 
the region of interest (pins a generons bnffer zone) at any time dnring the conrse of the 
simnlation are tagged. The entire simnlation is then reeled back to the nniform grid stage 
and for each of the top grid particles that were tagged as having passed throngh the region 
of interest, a set of smaller particles is generated for the snbgrid calcnlation. This is done 
in the following manner. Each of the tagged top grid particles dehnes a cnbical box of 
length It, eqnal to the interparticle spacing in the top grid. Allowing the snbgrid to be a 
factor of / smaller than the top grid, a virtnal grid of snbgrid particles is then generated 
with a spacing of Ig = h/f and any virtnal particle in a box dehned by a top grid particle 
is connted as a snbgrid particle. At this point the snbgrid particles constitnte a nniform 
grid which is fnlly sampled inside the snbgrid volnme and only partially sampled ontside it. 
Initial conditions for both the top grid and the snbgrid are then generated by pertnrbing the 
top grid and snbgrid particles away from their respective nniform grids. (Seeds identical 
to the seeds nsed to generate the hrst set of top grid initial conditions are nsed so that 
the initial conditions for the top grid in the mnlti-grid calcnlation will be identical to the 
initial conditions for the top grid alone.) 

The fnll mnlti-grid calculation is then run with the two sets of initial conditions, the 
top grid and the subgrid being evolved forward in time simultaneously. The important 
points to note are: (1) there is no “back-reaction” from the subgrid to the top grid (i.e. the 
top grid simulation runs completely unaware of the subgrid simulation), (2) the potential 
in the subgrid is computed using both the small particles in the subgrid and the potential 
from the top grid (i.e. the force held from the top grid acts as an external held on the 
subgrid simulation), and (3) unlike the top grid, the subgrid utilizes isolated boundary 
conditions so that subgrid particles may enter and exit the subgrid region over the course 
of the simulation. Additionally, the subgrids may either be kept stationary throughout the 
course of the simulation or they may be allowed to move (eg. to follow the growth of an 
object which has a large streaming velocity). 

An HPM simulation is constrained to use the same number of grid cells in both the 
top grid and subgrid simulations. Therefore by using a subgrid which is factor of / smaller 
than the top grid, the gain in spatial resolution in the subgrid region is necessarily a factor 
of /. The total number of particles used in the subgrid may, however, vary from that used 
in the top grid. Allowing an identical number of top grid and subgrid particles across a 
uniform grid, a subgrid which is a factor of / smaller than the top grid results in a mass 
per particle that is a factor of smaller than in the top grid. However, in a high density 
region of the simulation (eg. a cluster) it is oftentimes necessary to reduce the number of 
uniform grid particles in the subgrid compared to that of the top grid in order to remain 
within the available machine memory. 
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Very high resolution in an HPM simulation may be obtained by further nesting sub¬ 
grids within subgrids. Again, this is an iterative process. At the end of a two-level (top 
grid plus one subgrid) calculation, a region of interest is identihed within the subgrid. All 
of the subgrid particles which pass through that region of interest over the course of the 
simulation are tagged and the entire simulation (top grid plus subgrid) is reeled back to 
the uniform grid stage whereupon a second subgrid is generated within the hrst subgrid 
utilizing all the tagged particles from subgrid ^1. Using the same random number genera¬ 
tor seeds as were used previously, initial conditions for each of the grids (top grid plus two 
subgrids) are generated by perturbing each of the three sets of particles away from their 
respective uniform grids. Again, in the multi-grid calculation the grids are evolved forward 
in time simultaneously, there is no back-reaction from “parent” grid to “child” grid, and 
the force held from the “parent” grid acts as an external held for the force calculation in 
the “child” grid. 

Since the “child” particles experience high frequency power in their subgrid that their 
“parent” particles do not, the “child” particles will not move exactly in concert with 
their “parent” particles. However, the “child” particles do not stray very far from the 
general location of the “parent” since the extra high frequency power does not induce 
large streaming velocities. Should a “parent” particle exit/enter the region of a subgrid, 
it takes its “children” out of/into the region in a smooth manner. A simple consistency 
check involves comparing the number of “child” particles found inside a given subgrid at 
a particular time to the number of “parent” particles also within the subgrid at the same 
time. The ratio of these numbers should be of order the cube root of the ratio of the 
particle masses in the two grids (but will not be exactly equal to this value owing to the 
smooth manner in which the “child” particles enter and exit the subgrids). This is, indeed, 
the case and for the simulations presented here the ratio of the number of “child” particles 
to “parent” particles differs from the cube root of the ratio of the particle masses by an 
average of about 6%. 

Visual comparisons between the structures in “parent” and “child” grids show excel¬ 
lent agreement (see Fig. 1). However, due to the higher force resolution in the “child” 
grid, structures in the “child” grid tend to be more concentrated than the analogous, more 
poorly resolved structures in the “parent” grid. That is, the force in the “child” grid is 
not as soft as in the “parent” grid, allowing structures to collapse on smaller scales. 

3. THE SIMULATIONS 

Three multi-grid simulations of the formation of clusters in a standard CDM universe 
(Hq = 50 km/s/Mpc, Oq = 1, A = 0) were run. All simulations consisted of 3-level 
calculations: a top grid of length Ltop = 200h~^ Mpc, inside which was nested a subgrid 
of length Lsuhi = 33.3/i“^ Mpc, inside which was nested a subgrid of length Lsub 2 = 
8.3h~^ Mpc (comoving lengths). In all cases 256^ grid cells were used, resulting in a grid 
cell length of 32.6h~^ kpc (comoving) in the smallest, highest resolution subgrid. A total 
of 128^ particles were used in the top grid calculation, resulting in a mass per particle in 
the top grid of 1.06 x Owing to machine memory limitations, the uniform 

subgrids were constrained to fewer particles than the top grid. The particle mass in the 
low-resolution subgrids (subgrid ^1) was 3.90 x 10^^/i“^M© while in the high-resolution 


5 



subgrids (subgrid ^2) it was 4.88 x 10^At the end of the simulation there were 
8.3 X 10^ particles inside the high-resolution subgrid containing cluster 1, corresponding 
to a density of 179 particles per cubic megaparsec. For cluster 2, there were 7.3 x 10^ 
particles inside its high-resolution subgrid at the end of the simulation, corresponding to a 
density of 159 particles per cubic megaparsec. For cluster 3, there were 9.8 x 10^ particles 
inside its high-resolution subgrid at the end of the simulation, corresponding to a density 
of 211 particles per cubic megaparsec. 

The locations of the subgrids were specihed as follows. To begin, a top grid simulation 
was evolved from erg = 0.033 to ag = 1.0, where 


erg = 


— (8h ^Mpe) 
, P 


( 1 ) 


Identifying the hnal timestep of the top grid simulation as the present epoch (redshift, 
z, of 0), the simulation began at z = 29. This is a model which is somewhat under- 
normalized compared to the COBE observations (eg. Bunn & White 1997) and over- 
normalized compared to the abundance of rich clusters (eg. Bahcall & Cen 1993; White, 
Efstathiou & Frenk 1993; Eke, Cole & Frenk 1996; Viana & Fiddle 1996). To determine 
the present-day locations of rich clusters, the mass density held of the top grid at ug = 1.0 
was smoothed with a Gaussian hlter of length 1.5h~^ Mpe and the locations of peaks 
determined. From this smoothed density held the locations of three of the largest density 
peaks were selected as the centers of the hrst subgrids. For each of these subgrids, the top 
grid particles that passed through the subgrid region (plus 20% buher zones) were tagged, 
the simulations were reeled back to the uniform grid stage, initial conditions for 2-level 
calculations were generated, and the 2-grid simulations were then evolved from ug = 0.033 
to CTg = 1.0. The second, highest resolution subgrids were then chosen to be centered on 
the centers of mass of the clusters that formed in each of the hrst subgrids. Again, using 
the timesteps of the 2-level calculations, the particles in the hrst subgrids that passed 
through the regions of the second subgrids (plus buher zones) were tagged, the simulations 
were reeled back to the uniform grid stage, initial conditions for 3-level calculations were 
generated, and the 3-grid simulations were then evolved from ug = 0.033 to ug = 1.0. 

The clusters investigated here correspond to the very largest (i.e. most massive) of 
these objects that would typically form in CDM universes. They do not, therefore, repre¬ 
sent an unbiased, “average” sample of clusters but may correspond fairly with the “richest” 
clusters that would form in such universes. Tables 1,2, and 3 contain summaries of various 
properties of the clusters obtained from analyses of the highest resolution subgrids. At 
the end of the simulation all three clusters have masses of order 2 x 10 ^^/i“^Mq within 
the Abell radius (1.5h~^ Mpe), corresponding to of order 4 x 10^ particles in the highest 
resolution subgrids. 

Shown in Fig. 1 are grey-scale pictures of the clusters at ug = 1.0. The level of grey 
indicates the logarithm of the mass density along the line of sight in the projection and 
each projection has dimensions 8.3/i“^Mpc x 8.3/i“^Mpc x 8.3/i“^Mpc. That is, shown in 
Fig. 1 is a 2-dimensional compression of a 3-dimensional volume corresponding to the full 
volume of the highest resolution subgrid and each projection is centered on the center of 
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mass of the cluster. The top panels show the clusters at highest resolution (subgrid ^2) 
and the center panels show the clusters at lower resolution (subgrid ^1). The sizes of 
the pixels in the hgure correspond to the sizes of grid cells in the different subgrids and 
reflect their relative levels of resolution. In each case the high mass density in the inner 
regions of the clusters results in a “burned out” image, but in the outer regions of the 
clusters it is clear that there are many smaller galaxy-sized mass concentrations. (There 
are also smaller galaxy-sized mass concentrations in the inner regions of the cluster which 
are not visible in Fig. 1 due to the level of contrast. See §4.6.) The bottom panels in 
Fig. 1 show a comparison of the mass density along the line of sight in the high- and 
low-resolution subgrids. Specihcally, the grey-scale indicates the logarithm of the ratio of 
the mass density in subgrid ^2 to the mass density in subgrid ^1. The comparison is done 
at the same (low) resolution as subgrid ^1. Overall the comparison is excellent (i.e. the 
image is fairly flat at a moderate level of grey, indicating a density ratio of order unity). 
The largest discrepancy between the densities in the two subgrids occurs near the “edges” 
of the clusters where the density in subgrid ^2 is less than in the corresponding regions of 
subgrid ^1 (i.e. white pixels). The discrepancy is caused by the relative levels of numerical 
softening in the two subgrids, the force being softer on a larger scale in subgrid ^1 than 
in subgrid ^2, resulting in less concentrated structure in subgrid ^1 compared to subgrid 
^2. That is, in the outer regions of the cluster the high-resolution version is somewhat less 
dense than the low-resolution version since, on the whole, the cluster is more condensed in 
subgrid ^2 than it is in subgrid ^1. Spherically-averaged density prohles of the clusters 
computed using the two different subgrids are, however, in excellent agreement at large 
radii (see Fig. 5). 

4. RESULTS 

4.1 Mass Accretion 

The growth of each of the clusters was investigated through: [1] the mean infall 
distance of particles into the cluster, [2] the time evolution of the total mass of the cluster 
contained within the Abell radius and [3] the rate at which mass was accreted within 
the Abell radius as a function of time. In order to calculate each of these quantities 
the center of mass of each cluster is required. This was determined using the following 
iterative procedure. Starting with an initial center of mass given by the location of the 
corresponding peak in the smoothed mass density held of the top grid, all subgrid particles 
within a radius of 3.0/i“^Mpc were selected and the center of mass of those particles was 
computed. From this center of mass a new sphere of particles of radius 3.0/i“^Mpc was 
selected and a new center of mass computed. The process was repeated within a given 
subgrid until convergence was reached (of order 6 iterations). Note that the centers of mass 
for each cluster computed independently from the corresponding high- and low-resolution 
subgrids are identical. Also, over the course of the simulations the centers of mass of the 
clusters have low streaming velocities and they move total distances which are less than 
the mesh resolution of the top grid. 

Fig. 2 summarizes both the mean and maximum infall distances of particles over the 
course of the simulations. The points with error bars in Fig. 2 indicate the mean initial 
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distance of subgrid ^1 particles from the centers of mass of their respective clusters as a 
function of their distance from the cluster centers of mass at the end of the simulation. 
That is, the points indicate the mean streaming distance of particles present in the clusters 
at the end of the simulation as function of their distance from the center of mass at the 
hnal timestep. The error bars show one standard deviation. The open squares without 
error bars indicate the maximum initial distance of any one particle from the center of 
mass of its cluster versus its distance from that cluster at the end of the simulation. From 
this hgure, then, over the course of the simulations the mean distance traveled by particles 
found within the Abell radii at the end of the simulations is of order 12/i“^Mpc. This is 
not at all surprizing since the mass of the clusters within the Abell radius at ag = 1.0 is of 
order 2 x equal to the mass contained within a uniform critical density sphere 

of radius 12/i“^Mpc. The maximum infall distance, however, is of order 18/i“^Mpc. Thus, 
in order to follow all of the infall of mass into the clusters over the course of the entire 
simulation, a simulation that utilizes either constrained initial conditions or the excision of 
peaks from a large-scale structure simulation would require a volume of ~ 3 X lO^/i ^Mpc^ 
to be simulated at comparably high resolution. This is a factor of order 50 larger than the 
requisite volume for the highest resolution subgrid in the HPM calculation. 

The masses of the clusters contained within the Abell radius, M(r = 1.5/i“^Mpc), and 
the rates at which mass was accreted within the Abell radius, M(r = 1.5/i“^Mpc), were 
computed as a function of lookback time for each of the three clusters. The high-resolution 
subgrid particles were used for these calculations and the lookback time was computed by 
taking ag = 1.0 to be the present epoch. Results for the evolution of the total amount of 
mass within the Abell radius, normalized by the present-day mass of the cluster within the 
same (comoving) distance, are shown in Fig. 3. It is clear from this hgure that the details 
of the evolution of the clusters are somewhat different in each case but that all three gained 
of order 50% of their present-day mass within the past 5 to 8 Gyr . The details of the rate 
at which mass was accreted within the Abell radius are shown in Fig. 4. Here the rate at 
which mass was accreted by each cluster, M, is shown as a function of lookback time and 
redshift, normalized by fhe average rate at which the cluster accreted mass between 

z = 2 and z = 0. From this hgure it is clear that, although all three clusters have similar 
masses at the present, no single pattern of mass accretion dominates in the formation of 
the clusters. Cluster 1 shows a monotonic increase in mass accretion rate from z — 2 to 
z — 0.2, after which it accretes virtually no mass. Cluster 2, however, shows a monotonic 
decrease in the mass accretion rate from ^ = 2 to the present and cluster 3 forms via a 
mass accretion rate which is roughly constant. 

4.3 Density Profiles 

Spherically-averaged differential density prohles, p(r), are shown in Fig. 5 for each of 
the clusters at the end of the simulations. Squares indicate the density prohles obtained us¬ 
ing the high-resolution subgrid particles and triangles indicate the density prohles obtained 
using the low-resolution subgrid particles. Due to the numerical softening of the force on 
small scales the density is computed only on scales larger than two grid cells (65/i“^ kpc 
in the high-resolution subgrids and 260h~^ kpc in the low-resolution subgrids). Over the 
length scales for which the density can be computed in both the high- and low-resolution 



subgrids there is excellent agreement between the two calculations and the small-scale 
density prohle computed from the high-resolution subgrid is clearly a smooth continuation 
of the larger-scale density prohle computed from the low-resolution subgrid. Also, but for 
a suggestion of a battening in the density prohle of cluster 1, there is no turnover in the 
density prohles at small radii, which is as expected in purely dissipationless simulations. 
The apparent battening in the density prohle of cluster 1 may be partially numerical in 
origin as it occurs on scales corresponding to 3 grid cells and less. A higher resolution sim¬ 
ulation (eg. a third subgrid) would be required to determine whether the trend is indeed 
real or purely an artifact. A comparison of these density prohles and corresponding density 
prohles obtained from the particles in the top grid calculation is not warranted owing to 
the extremely poor resolution of the clusters in the top grid (of order 1000 particles in 
total and a force softer than Newtonian on scales smaller than the Abell radius). 

It is clear from Fig. 5 that the density prohles of the clusters are not well-ht by a single 
power law over all scales. This is to be expected since it has been shown previously that 
CDM halos are well-described by density prohles in which the logarithmic slope varies 
gently (eg. Dubinski & Carlberg 1991; Navarro, Frenk & White 1995, 1996ab; Cole & 
Lacey 1996; and Tormen, Bouchet & White 1997). On scales r < r 2 oo, where r 2 oo is the 
radius inside which the mean interior overdensity is 200, a good two-parameter ht to the 
density prohle is given by 


p{r) ^ 5c 

Pc x{l -t- x)‘^ 


( 2 ) 


(eg. Navarro, Frenk, & White 1996b) where pc is the critical density for closure of the 
universe, x = r/r^, and is a scale radius. Here 5c is a dimensionless characteristic 
density. By dehning the “concentration” of a halo to be c = r 2 oo/’'s, the two-parameter 
ht above can be reduced to a one-parameter ht through 


200 

3 [ln(l-F c) - c/(l-t-c)] 


( 3 ) 


(eg. Navarro, Frenk & White 1996b). 

Using the particles in the highest resolution subgrids, the values of r 2 oo (the “virial 
radius”) for the clusters were determined at erg = 0.67,0.83,1.0. The virial radii evolve 
relatively little from ug = 0.67 to erg = 1.0 and are of order 2h~^ Mpe for each of the 
clusters (see Tables 1, 2, and 3 for specihe values). Again using the particles in the highest 
resolution subgrids, the variation of the cluster overdensities with radius were evaluated 
on scales less than r 2 oo and results are shown by the points in Fig. 6. The solid lines in 
this hgure illustrate the best-htting density prohles of the form of equation (3) above. The 
values of the corresponding scale radii, Ts, are given in each of the panels of the hgure. But 
for a slight downturn in the small-scale density prohles of cluster 2 at ug = 0.67 and cluster 
3 at ug = 0.83, there is very good agreement between the simulated clusters and equation 
(3). Again, it is possible that the small-scale downturn is numerical in origin. The scaled 
density prohles in Fig. 6 are all fairly similar and in the case of clusters 1 and 2 the scale 
radius of the best-htting prohle evolves little from ug = 0.67 to erg = 1.0; however, for 
cluster 3 the value of changes appreciably (by a factor of order 2) over the same time 
period. 
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4.4 3-d Shapes S>i 2-d Ellipticity Distributions 

The evolution of the 3-dimensional and projected shapes of the clusters were computed 
at (jg = 0.67, 0.83, and 1.0 using the particles in the highest resolution subgrids. Since 
there is no hard “edge” to the clusters in terms of distinguishing those particles which are 
inside the cluster and those which are not, we shall define the boundaries of the clusters 
to be the virial radii (r 2 oo) for the following analyses. 

Using all particles within r 2 oo of the cluster centers of mass, the 3-dimensional cluster 
shapes were determined from a standard moment of inertia analysis that yielded the axis 
ratios h/a and c/a for each of the clusters (we define a > b > c). From the axial ratios a 
triaxiality parameter was computed for each of the clusters: 


T = 



(4) 


where T = 0 indicates a purely oblate object and T = 1 indicates a purely prolate object. 
We shall refer to objects with 0<T<l/3as being nearly oblate, those with 2/3 < T < 1 
as nearly prolate, and those with 1/3 < T < 2/3 as triaxial. Values of the cluster axial 
ratios and triaxiality parameters are listed in Tables 1, 2, and 3 for ug between 0.67 and 
1.0. From these tables it is clear that the evolution of the shapes of the clusters are quite 
different in each case. Although cluster 1 and cluster 2 are both nearly oblate at the end 
of the simulation, cluster 1 evolves from being triaxial at ug = 0.67 to being nearly oblate 
at CTg = 1.0 whereas the shape of cluster 2 changes little over the same period of time and, 
as a result, is always nearly oblate. Cluster 3, on the other hand, is nearly prolate at the 
end of the simulation but was nearly oblate at ug = 0.67. 

A more useful quantity for comparison of the evolution of the shapes of simulated 
clusters to observed clusters is the ellipticity projected on the plane of the sky. In the case 
of the simulations the projected ellipticity of the mass is the only quantity which can be 
computed reliably (i.e. without having to resort to assumptions about the degree to which 
mass would trace light). Given recent advances in gravitational lensing analyses of observed 
clusters, however, this seems a reasonable quantity to compute. That is, from analyses of 
the coherent weak distortion of the shapes of background galaxies due to an intervening 
gravitational lens cluster it is possible to constrain the ellipticity of the projected mass 
of clusters and, additionally, it is becoming clear that the smoothed light distribution of 
clusters traces the mass quite well (eg. Bonnet, Mellier & Fort 1994; Fahlman et al. 1994; 
Small et al. 1995; Kneib et al. 1996; Seitz et al. 1996; Squires et al. 1996ab; Small et al. 
1997). 

The ellipticities of the clusters as projected on the sky, e = 1 — 6/a, were computed 
using all particles in the the highest resolution subgrids that were located within a distance 
of r 2 oo of the cluster centers of mass. The probability, T'(e), of observing a given projected 
ellipticity for a given cluster was computed by viewing each cluster from 500 random 
orientations and assembling an appropriately normalized probability distribution function. 
Results are shown in Fig. 7 for all 3 clusters at ug = 0.67, 0.83, and 1.0. Again, as with the 
evolution of the triaxiality parameter, each cluster exhibits its own particular evolution 
in projected shape. Cluster 1 evolves toward being, on average, significantly flatter in 
projection at ug = 1.0 than it was at ug = 0.67. Cluster 2, on the other hand, remains 
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approximately the same projected shape over the same time period and cluster 3 evolves 
toward being signihcantly rounder in projection on average. The median ellipticities for 
each of the clusters, Cmed, are listed in Tables 1, 2, and 3. 

4.5 Substructure 

The redshift dependence of the fraction of observed clusters having a signihcant 
amount of substructure is potentially a good indicator of the value of density parame¬ 
ter. This is due to the fact that in a universe in which Oq < 1 density fluctuations cease 
to grow at redshifts of order — 1 and, so, moderate to low redshift clusters in critical 
density universes are expected to contain substantial amounts of substructure on average 
while in low density universes the clusters should be much more regular. Wilson, Cole & 
Frenk (1996) have explored the possibility of using observations of weak lensing by clusters 
to discriminate between universes with low and critical values of Oq via a quantihcation of 
cluster substructure from the weak shear held. Although initially optimistic, the situation 
has become more murky recently with the realization that some of the simulated clusters 
used in the analysis were inadequate. 

The simulations discussed here are restricted to a critical density universe and, hence, 
we do not investigate the explicit Oq dependence of cluster substructure with cosmological 
epoch (this analysis will be performed in future simulations). Rather, we have investigated 
the evolution of substructure in the clusters from erg = 0.67 to ag = 1.0 in order to asses 
the degree to which substructure is erased over this period of time. 

There are numerous methods by which cluster substructure can be quantihed but 
here we restrict the analysis to the Dressier-Shectman A statistic (Dressier & Shectman 
1988). This choice is made based on the results of Pinkey et al. (1996) who have subjected 
many substructure tests to thorough analysis and conclude that by and large the Dressler- 
Shectman test tends to be the most sensitive to substructure. 

The A statistic is dehned by 



i=l 


sf 



vf + (cTi 



(5) 


where N is the number of galaxies in the cluster and Vi and are, respectively, the mean 
velocity and the velocity dispersion of the Woe nearest neighbors to each galaxy. The 
sensitivity of the A statistic is dependent upon the number of neighboring galaxies used 
in the analysis and Bird (1995) hnds the test to be most sensitive for Woe = VN- The 
statistic is a measure of the correlation between the (projected) locations of the galaxies in 
the cluster and their velocities. In the case of uncorrelated positions and velocities A ~ 1. 
The quantihcation of substructure for a given cluster using only the computed value of A 
is insufficient, however, and in order to assess the likelihood of real substructure within 
the cluster Monte Carlo simulations must be performed. Additionally, Crone, Evrard & 
Richstone (1996) have pointed out that since the value of A is not independent of N, 
the total number of galaxies used in the analysis, in order to compare either observed or 
theoretical clusters to one another it is necessary in the analyses to select an identical 
number of galaxies for each cluster. 
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The simulations presented here are of insufficient resolution to resolve the dark matter 
halos of individual galaxies (see §4.6 below) and we instead investigate the substructure in 
the mass distribution. (For all intents and purposes this is the nature of the substructure 
identihed via weak lensing, though there is clearly good correspondence between “lumpi¬ 
ness” in the mass and galaxy distributions in lensing clusters.) Substructure in the clusters 
was evaluated in the following way. Since the specihc value of A for a given cluster will 
depend on the angle from which it is viewed in projection, each cluster was viewed from 500 
random directions. For each viewing angle, the projected locations of a randomly selected 
subset of the particles was used to compute a value of A. Additionally, for each viewing 
angle a Monte Carlo value of the statistic (Arand) was computed by randomly shuffling the 
velocities of the particles amongst their positions. Mean values of A and Aj-and, along with 
their formal standard deviations, were then computed from the 500 individual values. For 
each cluster a total of A = 1024 particles were randomly selected from all particles within 
the virial radius (r 2 oo)- Different random sets of particles were used for each viewing angle 
and for the timesteps corresponding to as = 0.67, 0.83, and 1.0. Following Bird (1995), 
Aloe was taken to be \/N = 32. 

Tables 1, 2, and 3 list the mean values of the A statistic for each cluster computed 
from the 500 random viewing angles, the mean values computed from the 500 Monte 
Carlo position-velocity shuffles (A^and), and the corresponding 1-a errors. Also listed are 
the ratios A/Arand, which indicate that all three clusters contain signiheant amounts of 
substructure at each of the three epochs, as- That is, within the virial radius substructure 
in the mass distribution of the clusters is not completely erased by the end of the simulation. 
In the case of cluster 2, the degree of substructure over and above the expectations of 
random is roughly constant at a ~ 3-a level from as = 0.67 to as = 1.0. In the case of 
clusters 1 and 3, there is signihcantly less substructure at as = 1.0 than at as = 0.67. 
However, the “erasure” of substructure in these two clusters over this time period is not 
quite monotonic. Cluster 1 is a bit “lumpier” at as = 0.83 than it is at either as = 0.67 
or Us = 1.0, while cluster 2 is a bit smoother and less “lumpy” at as = 0.83 than it is at 
either as = 0.67 or as = 1.0. 

Bubble plots of the A test are very helpful for illustrating visually both the location 
and amount substructure in a cluster. Fig. 8 shows bubble plots for each of the clusters 
at the end of the simulation. The viewing angle for the projection of each cluster was 
chosen to be an angle for which the specihc value of A was identical to the mean value 
of the 500 random orientations. The dots in Fig. 8 show the spatial locations of the 1024 
randomly selected mass points used in the A statistic analysis. The circles, all of which are 
centered on dots, have been scaled to have radii proportional to Si (see equation 5 above). 
The larger the circle, the larger is the local deviation of the mean velocity and/or velocity 
dispersion from the global cluster value. For clarity, circles are drawn around only those 
mass points for which Si > 2A (i.e. regions of most signiheant substructure). Far from 
being smooth blobs of mass, the clusters are all clearly “lumpy”, each having of order 4 or 
5 signiheant sub-lumps within the virial radius. 

4.6 Galaxy Halos 

It is often thought, erroneously, that purely dissipationless simulations are inadequate 
to study the dark matter halos of galaxies in a simulated cluster environment because all 
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small dark matter concentrations are destroyed by pnrely nnmerical effects as they orbit 
throngh the clnster. However, Bromley et al. (1995) have demonstrated that provided the 
force resolntion is high galaxy-sized dark matter halos will snrvive many orbits throngh 
the potential of a large clnster and are not destroyed by pnrely nnmerical effects. 

In the case of this work, the stndy of individnal galaxy-sized dark halos in the clnsters 
is not a reasonable goal since the length scale over which the force is non-Newtonian is 
somewhat too large, even in the case of the highest resolntion snbgrids. However, it can 
be seen in Fig. 1 that there are nnmerons galaxy-sized concentrations of dark matter in 
the onter regions of the clnsters and the clnsters are not merely smooth blobs of mass. 
(The Inmpiness of the mass distribntion is borne ont in part by the snbstrnctnre analysis 
above.) Dne to the high particle density in the central regions and the choice of contrast 
level, existing galaxy-sized concentrations in the inner regions of the clnsters are not visible 
in Fig. 1. 

Althongh it is clearly an inadeqnate method for the generation of a highly accnrate 
catalog of galaxy-sized dark halos, a simple friends-of-friends algorithm was nsed to gen¬ 
erate catalogs of gronps of particles in the highest resolntion snbgrids. The gronps were 
selected to have overdensities d ^ 1000 (typical of the overdensity of the Inminons region 
of a bright galaxy) and masses ^ 10^^h~^MQ (20 particles or more). At the end of the 
simnlations, ~ 300 snch objects were fonnd within the Abell radii of clnsters 1 and 3, and 
~ 200 were fonnd within the Abell radins of clnster 2. However, owing to the large scale 
over which the force is softer than Newtonian, the central 0.5/i“^Mpc of each clnster is 
dominated by a single hnge “halo” of mass 

By nesting yet another snbgrid within the highest resolntion snbgrid (i.e. by perform¬ 
ing a 4-level calcnlation), the effective length scale over which the force is softer than 
Newtonian will be rednced significantly. Within the Abell radins of a clnster it will then 
be possible to resolve confidently gronps of particles that may be fairly associated with 
the dark matter halos of individnal galaxies and to eliminate the artificial overmerging of 
halos in the central region. Snch analyses will be performed in fntnre simnlations. 

5. DISCUSSION 

Using a hybrid N-body code in which high mass and spatial resolntion can be obtained 
in small regions of a very large total simnlation volnme, the formation of three massive 
clnsters was investigated. The clnsters were chosen to be typical of the most massive 
clnsters that wonld be present in a standard CDM nniverse at an epoch corresponding to 
CTg = 1.0. At highest resolntion, the clnsters consisted of ~ 4 x 10^ particles within the 
Abell radins at as = 1.0. Althongh the clnsters share similar properties at the end of the 
simnlation, the details of their formation histories are qnite different. 

The properties which the clnsters share are: 

- formation within the same large compntational volnme (200^/i“^Mpc^) 

- masses of M(r < 1.5/i“^Mpc) ~ 2 x 10^^h~^MQ at the end of the simnlation 

- similar spherically-averaged density profiles which are well-fit by eqnation (2) 

- similar valnes of the virial radins, r 2 oo, and similar valnes of the scale radins, r^, at 
the end of the simnlation 
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- accretion of ~ 50% of the mass present at the end of the simnlation within the past 

5 to 8 Gyr 

In terms of the details of the formation history of the clnsters, however, each clnster 
exhibits markedly individnal behavior and no single pattern dominates in the evolntion of 
the following clnster properties: 

- mass accretion rate, M 

- scale radins, 

- three-dimensional shape (triaxiality parameter, T) 

- two-dimensional (projected) ellipticity 

- snbstrnctnre (Dressler-Shectman A statistic) 

Based on the very small nnmber of clnsters presented here it is difhcnlt to make 
any statistically-sonnd conclnsions abont clnster formation and evolntion. However, the 
problems with standard CDM notwithstanding, there are certainly some interesting things 
to be noted. The nnmerical clnsters are extremely massive and, so, correspond to rich 
clnsters. Rich clnsters, being the brightest and most massive, are likely to be those which 
will be stndied observationally over the widest range of redshifts and are the objects from 
which it is hoped that cosmological constraints will arise. Given the markedly different 
formation histories of the clnsters, a qnestion raised by this investigation is the degree 
to which observations of a sample of rich clnsters covering a wide range of redshifts can 
provide stringent cosmological constraints. 

Althongh the nnmerical clnsters stndied here do not constitnte a statistically large 
sample, it is clear that snch large clnsters do not form a simple 1-parameter family as far 
as their evolntion history is concerned, even thongh the clnsters have similar masses and 
spherically-averaged density prohles at the present day. The nse of observations of clnster 
evolntion to constrain cosmological models may yet be viable bnt the cantions snggestion 
from this work is that this may not be completely straightforward. Gonsiderably more work 
on the details of the formation history of clnsters in varions models of strnctnre formation 
is necessary in order to determine both the degree to which cosmological conclnsions can 
be drawn from observations of clnsters at different epochs and also the reqnisite size of 
a sample of observed clnsters which wonld insnre those cosmological conclnsions to be 
statistically reliable. 
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FIGURE CAPTIONS 

Fig. la: Grey-scale images of clnster 1 at ag = 1.0. The top panel shows the logarithm of 
the mass density along the line of sight in the high-resolntion snbgrid (snbgrid ^2) and 
the center panel shows the same for the low resolntion snbgrid (snbgrid ^1). The pixel 
sizes in the image are eqnal to the grid cell sizes in the snbgrids and the projection is of a 
8.3/i“^Mpc X 8.3/i“^Mpc x 8.3/i“^Mpc cnbe, centered on the center of mass of the clnster. 
The bottom panel shows a comparison of the line of sight mass density in the two snbgrids, 
where the grey scale indicates the logarithm of the mass density in snbgrid ^2 divided by 
the mass density in snbgrid ^1, compnted at the (low) resolntion of snbgrid ^1. Overall, 
the mass density in the two snbgrids compares well (see text). 

Fig. lb: Same as Fig. la, bnt for clnster 2. 

Fig. Ic: Same as Fig. la, bnt for clnster 3. 

Fig. 2: The infall distance of particles into the clnsters. Solid points with error bars show 
the mean distance streamed since the beginning of the simnlation as a fnnction of the 
hnal location of the particles, compnted relative to the clnster centers of mass. Frror bars 
indicate one standard deviation. The open sqnares show the maxim um initial distance of 
any one particle from the center of mass of its clnster as a fnnction of its location at the 
the end of the simnlation. 

Fig. 3: Mass of the clnsters contained within the Abell radins as a fnnction of lookback 
time. The contained mass has been normalized by the contained mass at the end of the 
simnlation. All three clnsters have masses ~ 2 x 10^^/i“^Mpc within the Abell radins at 
the end of the simnlation (see Tables 1, 2, and 3). 

Fig. 4: Rate at which mass is accreted, M, within the Abell radins as a f un ction of lookback 
time. The accretion rate is normalized by the mean rate, which mass is accreted 

within the Abell radins between lookback times corresponding to x = 2 and z = 0. 

Fig. 5: Spherically averaged density prohles for each of the clnsters, evalnated at ag = 1.0. 
Open sqnares indicate p(r) compnted nsing the high-resolntion snbgrid particles; hlled 
triangles indicate p(r) compnted nsing the low-resolntion snbgrid particles. 
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Fig. 6: Cluster overdensities as a function of radius (scaled by r 2007 the virial radius) for 
CTg = 0.67, 0.83 and 1.0. Results for cluster 1 are shown in the three panels on the left, 
results for cluster 2 are shown in the three central panels, and results for cluster 3 are 
shown in the three panels on the right. The solid line indicates the best-ht density prohle 
of the form of equation (2), with the corresponding value of the scale radius, rg, indicated 
in each of the individual panels of the hgure. 

Fig. 7: Projected ellipticity distributions for each of the clusters as a function of erg. The 
dashed line indicates erg = 0.67, the dotted line indicates ug = 0.83, and the solid line 
indicates ug = 1.0. 

Fig. 8: Bubble plots of the A test for each cluster at ug = 1.0. The linear scale of the 
projection is 5/i“^Mpc by 5/i“^Mpc. Dots indicate the spatial location of the mass points 
used in the evaluation of A and the circles have radii proportional to 5i (see text). The 
degree of substructure apparent in these projections is indicative of the mean value, based 
on the results of 500 random viewing angles. 
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Table 1: Properties of cluster 1 (results of highest resolution subgrid) 



(Tg — 0.67 

(Ts = 0.83 

O 

II 

GO 

b 

M{r < 1.5/?, ^Mpc) 

8.9 X 10’^/?,-’AfQ 

1.8 X 

1.9 X 101®/?,-iMq 

Nj,(r < 1.5/?,^^Mpc) 

1.8 X 10® 

3.7 X 10® 

3.9 X 10® 

^200 

1.60/?“^Mpc 

2.18/?,-^Mpc 

2.27/?,-iMpc 

h/a 

0.87 

0.91 

0.97 

c/a 

0.76 

0.76 

0.68 

T 

0.60 

0.38 

0.12 

^med 

0.30 

0.27 

0.38 

A 

1.75 ± 0.12 

1.86 ±0.12 

1.30 ±0.09 

^rand 

1.05 ± 0.07 

0.97 ±0.07 

1.04 ±0.07 

^rand 

1.67±0.16 

1.92 ±0.18 

1.25 ±0.12 


Table 2: Properties of cluster 2 (results of highest resolution subgrid) 



(Tg — 0.67 

(jg — 0.83 

O 

II 

GO 

b 

M{r < 1.5/?, ^Mpc) 

1.2 X 

1.5 X 

1.7 X 101®/?,-iMq 

Nj,{r < 1.5/?,-iMpc) 

2.4 X 10® 

3.0 X 10® 

3.5 X 10® 

^200 

1.83/?“^Mpc 

2.05/?-iMpc 

2.08/?,-iMpc 

h/a 

0.98 

0.98 

0.97 

c/a 

0.76 

0.76 

0.79 

T 

0.10 

0.08 

0.15 

^med 

0.30 

0.30 

0.26 

A 

1.45 ± 0.10 

1.43 ±0.10 

1.49 ±0.10 

^rand 

1.06 ± 0.08 

1.06 ±0.07 

1.07 ±0.07 

^rand 

1.38 ± 0.14 

1.35 ±0.13 

1.39 ±0.13 


Table 3: Properties of cluster 3 (results of highest resolution subgrid) 



(Tg — 0.67 

(Tg = 0.83 

O 

II 

GO 

b 

M{r < 1.5/?-iMpc) 

1.2 X 10’®/?,-’AfQ 

1.5 X 

2.1 X 101®/?,-iMq 

Nj,{r < 1.5/?,“^Mpc) 

2.4 X 10® 

3.0 X 10® 

4.3 X 10® 

*’200 

1.87/?,-iMpc 

2.06/?,“^Mpc 

2.24/?,-iMpc 

h/a 

0.96 

0.94 

0.87 

c/a 

0.70 

0.77 

0.86 

T 

0.16 

0.31 

0.93 

^med 

0.30 

0.30 

0.21 

A 

1.76 ± 0.12 

1.39 ±0.09 

1.48 ±0.09 

^rand 

1.02 ± 0.06 

1.08 ±0.07 

1.06 ±0.07 

^rand 

1.72 ± 0.15 

1.29 ±0.12 

1.40 ±0.12 
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Figure la 
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Figure lb 
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Figure Ic 
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